#!/usr/bin/R
	
#set model parameters
	t.horizon=1000 #simulation length (hrs)
	dt=.1 #timestep length (hrs)
	ET=0 #evapotranspiration rate (m/hr) ###consider functional relationship to temperature???###	
	PP=0 #precipitation rate (m/hr) ###forecasted???###

#intialize variables
	source('pondSim.module.R')
	source('pond.params.R')
	n.ponds=length(area)
	n.boxs=length(box.elev)
	n.times=as.integer(t.horizon/dt) # number of timesteps
	h=array(0,dim=c(n.ponds,n.times)) #water surface elev above MSL
	plant.inflow=array(0,dim=n.times) #inflow to first pond (cms)
	pond.inflow=array(0,dim=c(n.ponds,n.times)) #net inflow to each pond (cms)
	pond.outflow=array(0,dim=c(n.ponds,n.times)) #net outflow from each pond (cms)
	sea.level=0 #set constant sea surface elev
	q=array(0,dim=c(n.boxs,n.times)) #flow past each connection (cms)
	q.temp=array(0,dim=n.boxs) #dummy inflow for RK routine
	sea.level=0 #constant sea surface elev
	K=array(0,dim=c(3,n.ponds)) #dummy variables for RK routine
	h[1:n.ponds,1]=initial.elev #set initial conditions
	#colors##########################################
	library(fields,verbose=FALSE)
        my.colors=array(tim.colors(),dim=(n.times+1))
	#colors##########################################
	
#perform simulation
	#png(file='images/1.png',width=500,height=400,units='px',bg='transparent')
	barplot(h[,1],col=my.colors[1],ylim=c(0,5)) #plot initial conditions
	#dev.off()
	if (n.times==1) stop()
	for (t in 1:(n.times-1)) {
		#RK STEP 1
		foo=find.flow(h[,t], n.boxs, box.links, box.pump, box.elev, box.width, box.type, roughness)
			pond.inflow=foo[[1]]+(PP*area)
			pond.outflow=foo[[2]]+(ET*area)
			pond.inflow[1]=pond.inflow[1]+plant.inflow[t] #add plant inflow to first pond
			K[1,]=dt/area*(pond.inflow-pond.outflow)
			h.temp=h[,t]+(1/3)*K[1,]
		#RK STEP 2
		foo=find.flow(h.temp, n.boxs, box.links, box.pump, box.elev, box.width, box.type, roughness)
                        pond.inflow=foo[[1]]+(PP*area)
                        pond.outflow=foo[[2]]+(ET*area)
                        pond.inflow[1]=pond.inflow[1]+plant.inflow[t] #add plant inflow to first pond
                        K[2,]=dt/area*(pond.inflow-pond.outflow)
                        h.temp=h[,t]+(2/3)*K[2,]
		#RK STEP 3
		foo=find.flow(h.temp, n.boxs, box.links, box.pump, box.elev, box.width, box.type, roughness)
                        pond.inflow=foo[[1]]+(PP*area)
                        pond.outflow=foo[[2]]+(ET*area)
                        pond.inflow[1]=pond.inflow[1]+plant.inflow[t] #add plant inflow to first pond
                        K[3,]=dt/area*(pond.inflow-pond.outflow)
                        h.temp=h[,t]+(1/4)*K[1,]+(3/4)*K[3,]
		#RK STEP 4
		h[,t+1]=h.temp
		#png(file=paste('images/',t+1,'.png',setp=''),width=500,height=400,units='px',bg='transparent')
		barplot(h[,t+1],col=my.colors[t+1],ylim=c(0,5))
		#dev.off()
	}
